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We report results on the interrelation between driving force, roughness exponent, branching and 
crack speed in a finite element model. We show that for low applied loadings the crack speed 
reaches the values measured in the experiments, and the crack surface roughness is compatible with 
logarithmic scaling. At higher loadings, the crack speed increases, and the crack roughness exponent 
approaches the value measured at short length scales in experiments. In the case of high anisotropy, 
the crack speed is fully compatible with the values measured in experiments on anisotropic materials, 
and we are able to interpret explicitly the results in terms of the efficiency function introduced by us 
in our previous work [A. Parisi and R. C. Bah, Phys. Rev. B 66(16) 165432 (2002)]. The mechanism 
which leads to the decrease of crack speed and the appearence of the logarithmic scaling is attempted 
branching, whilst power law roughness develops when branches succeed in growing to macroscopic 
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I. INTRODUCTION 

There are two questions in the field of fracture me- 
chanics to which intense research has been devoted in 
recent years. The first concerns the terminal speed pre- 
dicted by the continuum theory^ which does not match 
the maximum speed measured in experiments Ter- 
minal crack speeds are usually found to vary in a range 
between about 90 % of the Rayleigh speed in anisotropic 
materials f^&i down to values as low as 33% for more 
isotropic materials The Rayleigh surface wave speed 
vr is the terminal crack speed expected in the contin- 
uum theory; beyond it the continuum solutions show 
compressive sign to what should be the crack opening 
stress component ahead of the crack tip.^^ The exper- 
imental maximum speed is usually accompanied by tip 
branching which is also not well explained. What con- 
trols the instability which leads to branching? What sets 
the terminal crack speed? Both the continuum theory 
and computer investigations have suggested that the an- 
swer to these questions lies in the mechanism through 
which energy is dissipated at the crack tip.^ 

The second question involves the roughness exponent 
governing the self-affine scaling of height fluctuations of 
fracture surfaces. Recent measurements^^ have shown 
how this quantity has a universal behaviour in that it 
takes a value of 0.5 at "short length scales" (in experi- 
ments this usually corresponds to nanometer scales) and 
a value of 0.75 at "large scales", the two regimes be- 
ing separated by a material dependent crossover length. 
In addition, a logarithmic scaling of fracture surfaces 
has been theoreticall}ii2ii^ and experimentally^^ found in 
the limit of quasi-static crack advance in brittle materi- 
als. Despite numerous attempts to describe the universal 
character of such behaviour, the question of what controls 
the value of the roughness exponent and why fractures 
tend to grow rough surfaces is still open. 

The two problems have repeatedly appeared to be con- 
nected. Experiments performed on polymethilmethacry- 



late (PMMA) have shown that beyond a critical crack 
speed, the crack tip starts to oscillate leading to the 
formation of structures on the crack surface and to a 
departure of the crack speed from that expected from 
the continuum theory. The same phenomenology was 
observed in simulations, where departure from steady 
state propagation for cracks exceeding a threshold speed 
was observed, with zig-zag motion and formation of 
microstructures.^^ Molecular dynamics simulations of 
crystalline silicon showed that cracks can dissipate large 
amounts of energy, up to seven times the energy needed 
to create a smooth surface as estimated in the framework 
of the continuum theory, the suggestion being that this 
energy goes into lattice oscillations^ 

The idea that the energy available does not all go into 
fracture work is not new. Analytical studies of planar 
crack advance in a lattice by Slepyan had already shown 
that the presence of the lattice leads to an important 
excess of energy being radiated from a crack at both low 
and high speeds^ and that crack propagation at low 
speeds is unstable^^. More recently, Marder and LiuiS 
have studied a lattice model for fractures, concluding that 
it is lattice oscillations which limit the range of possible 
crack speeds^2iS& 

In a recent paper^^ we have shown how such energy 
radiation at the crack tip due to phonon emission is the 
crucial mechanism for crack propagation. The intensity 
of the radiated energy is a function of the crack speed, 
ruling out stable crack growth at low speeds. The depen- 
dence of the crack speed is in agreement with the analyt- 
ical results of Slepyan. ^^'^^ The instability at low crack 
speeds also corresponds to what suggested in Ref. 18 and 
to the results obtained by Marder and Liui^^i^ Our re- 
sults suggest a way to relate the intensity of the radiated 
energy directly to the phonon band structure. The re- 
lation between this approach and the results of Marder 
et a/.^'^^'^^ is discussed in Ref. |^ The mechanism pro- 
posed has been investigated in the limiting case of planar 
cracks, but the analysis becomes more difficult in pres- 
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ence of branching. However there is clear simulation evi- 
dence in two dimensions that crack branching is sensitive 
to the lengthscale at which the continuum description 
breaks down.^^ 

The mechanism that leads to branching has not yet 
been fully unravelled. The first hint to the understand- 
ing of such phenomenon came from Yoffe^^ who showed 
that beyond a critical crack speed, the hoop stress has a 
maximum at a definite angle with respect to the direction 
of crack propagation. This has attracted considerable 
discussion as experiments have shown branching at crack 
speeds which differ from the prediction of Yoffe.^'^'^^ Sim- 
ulations and experiments however seem to agree that the 
mechanism of branching could be connected to the ter- 
minal crack speedA2SiSi 

In this work we will show that in our simulations the 
main mechanism that limits the crack speed for a free 
running crack is attemped branching. Cracks constrained 
on a plane reach high speeds compatible with the speeds 
measured in highly anisotropic materials, which we are 
able to explain in terms of the efficiency function intro- 
duced in our previous work.^^ We will also show that 
branching is responsible for the roughness exponent of 
the fracture surfaces: attempted branching roughens the 
crack surface with a logarithmic scaling first and, when 
macroscopic branches develop, with a roughness expo- 
nent close to the value measured at short length scales. 
At the same time, the crack speed is drastically reduced 
due to the attempted branching mechanism to values 
comparable with those measured in experiments, and 
gradually rejoins the efficiency description for high load- 
ings. 

The model used for these simulations has been exten- 
sively described in our previous work,^^ and its distinc- 
tive features are reviewed in section ^ In section IIIII 
we discuss the origins and effect of both disorder and 
anisotropy, and show how to implement disorder in the 
model and control the driving energy. Results on the 
crack speed for both planar and non-planar cracks are 
described in section II VI whilst the scaling of crack sur- 
face roughness for different driving regimes is discussed in 
section|3 An attempt to simulate disconnected fractures 
is reported in section IVTI and finally we draw conclusions 
in section IVTTI 



II. THE MODEL 

The model used in these simulations is an application 
of the finite element scheme and has been fully described 
in Ref. The elastodynamic description is obtained 
by discretizing space on an fee lattice and connecting 
neighbouring lattice points with (non-filling) tetrahedral 
elements. We then solve the Euler-Lagrange equations 
obtained by using the discretized form of the Lagrangian 
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FIG. 1: The efficiency for planar crack propagation in our 
fee lattice model. The overall increase of E{v) with the crack 
speed still remains to be quantitatively understood, but the 
drop of the efficiency at the Rayleigh speed as well as the 
fine structure are well understood in terms of vibrational 
resonances. 

of continuum elasticity: 

V t 

where is the displacement field at site v and at is 
the stress tensor at element t, the index v spanning all 
lattice points and the index t spanning all tetrahedral 
elements. The stress tensor is related to the displace- 
ment field by the standard tensorial relation of contin- 
uum elasticity a= ATr(Vu) 1 + /i[Vu + (Vu)^] and ft^ 
is a volume element related to the discretization scheme 
used. No mechanism of dissipation is active, but due to 
the discreteness introduced waves are radiated from the 
crack tip with an intensity which is non linear in the crack 
speed, revealing a selection rule for the crack speed itself. 

Rupture is simulated by irreversibly setting the elas- 
tic constants A and /i to zero in any tetrahedral element 
where the elastic energy stored within becomes greater 
than a pre-determined fracture energy. This is not the 
only possibility for a rupture criterion: quite recently, 
Heizler, Kessler and Levine^^'^^ have studied the conse- 
quence of having a set of continuous force laws as op- 
posed to the usual piecewise discontinuous force laws we 
use. Their results suggest that the nature of the force 
law could change the stability limit for high speed crack 
motion and, at least in mode-I fractures, for low crack 
motion. Using finite elements, a continuous force law 
was used by Johnson^^ to study the extension of the pro- 
cess region for moving cracks. Needleman implemented 
crack rupture in a finite element model using surface 
decohesion.^^ The model was used both in the study 
of crack growth in brittle solids^^ and in the study of 
interfacial crack growth^2i2i It was found that a crack 
constrained to a plane can approach the Rayleigh speed, 
whilst when unconstrained, crack branching and reduc- 
tion of the terminal crack speed were found. A similar 
behaviour is found in our simulations. 
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Our model has been used in our previous work in 
the simphfied case of two dimensions, to simulate pla- 
nar cracks running at fixed speed. In the present arti- 
cle, we present results of full three-dimensional simula- 
tions, in which both the condition of fixed crack speed 
and planarity are released, and the crack advance is only 
controlled by a Griffith criterion. In the continuum de- 
scription, the Griffith criterion states that a crack will 
advance only if the macroscopic energy delivered to the 
crack tip Gm('^7^) exceeds the fracture work 270 neces- 
sary to create new surface: 

GM{v,t) >27o 

If GM{v^t) > 270, the excess of energy is usually con- 
sidered the source of kinetic energy for crack advance, so 
if the loading is just sufficient to have the crack prop- 
agate, the crack should advance quasi-statically. Both 
simulations and experiments show that this is not the 
case: cracks do accelerate rapidly towards a limiting 
crack speed which depends on the material A^M^ 

The solution of the puzzle is the presence of an alter- 
native mechanism of energy dissipation due to the dis- 
creteness of matter at the microscopic level, a mechanism 
which is not included in the continuum elastodynamic de- 
scription. Due to the discreteness, the macroscopic en- 
ergy release rate Gm('^,^) can be equated to the sum of 
two microscopic contributions: 

Gm{v, t) = G\,r{v, t) + Gph(v, t) 

where G\^^{y^t) and G^h.{y^t) are respectively the break- 
age energy release rate, which is the portion of the avail- 
able energy going into fracture work, and the phonon 
energy release rate which is the portion of the available 
energy radiated as phonons. The breakage energy release 
rate can be expressed as: 

G^,r{v,t)=E{v)GM{v,t) (2.1) 

by introducing the efficiency E{v). In the case of a strip 
geometry with fixed displacements at the top and bottom 
boundaries, the macroscopic energy release rate is a con- 
stant independent of the crack speedi Gm('^,^) = G^. 
Hence, for this special case the efficiency becomes the sole 
source of velocity dependence, thus separating the effect 
of local discreteness from the effect of the macroscopic 
external loading. Moreover, the efficiency E{v) has been 
shown to depend only on the lattice geometry and crack 
speed, being local to the crack tip and independent of 
the macroscopic dynamical history.^^ 

The speed dependence of the efficiency for our fee lat- 
tice model is shown in fig. ^ Although the overall in- 
crease of E{v) with the crack speed still remains to be 
quantitatively understood, the drop of the efficiency at 
the Rayleigh speed as well as the fine structure are well 
understood and correspond to resonant emission, when 
emitted waves have a group velocity matching the crack 
speed itself. In particular, we expect the drop at the 



Rayleigh speed to be a feature common to all materials, 
as that resonance arises in the continuum limit. The effi- 
ciency E{v) for energy transfer into bond breakage is well 
below unity even for zero speed, as can clearly be seen in 
figure nj This is because when a lattice element (tetra- 
hedron) breaks at the crack tip, others around recoil dy- 
namically and that recoil energy is ultimately radiated 
as sound waves. 

A crack will only advance if the energy available at 
the crack tip is sufficient to create new surface. This 
translates the Griffith criterion in presence of discreteness 
into a condition of steady crack growth given by: 

Gbr(^,t) =270, (2.2) 

the difference from the macroscopic energy release rate 
being converted into phonons. In the case of a strip 
geometry with fixed displacements at the top and bot- 
tom boundaries, the threshold to initiate crack advance 
is given by G^ = 27o/£^(0) and therefore if the loading is 
maintained only cracks with speed such that E(v) = £^(0) 
can propagate, leading io v 0.88 1'i^ from the data of 
figure n We have further argued in Ref. i21| that only 
the more limited regions where dE/dv < should be 
sustainable. 

One can of course load a sample above the quasistatic 
threshold, with 

where e > 1. In this case the condition for sustained 
crack propagation becomes 

E{v) = E{0)/€ (2.3) 

and by reference to fig. we see that the higher loading 
can sustain higher crack speeds, as might be expected. 

As the efficiency only depends on the crack speed, by 
using eq. ^2.1^ we can calculate Ghr{v^t) for any macro- 
scopic loading for which Guiy^ t) is known, and use it in 
eq. (|2.2j) to get the corresponding allowed crack speeds. 

III. THE ROLE OF DISORDER AND 
ANISOTROPY 

Discretization naturally introduces preferred directions 
in space and therefore simulations are characterized by 
some level of anisotropy. In the absence of disorder, a 
square (cubic) lattice acts as a planar guide for the ad- 
vancing crack. This phenomenon, known as lattice trap- 
ping^ was predicted by Thomson, Hsieh and Rama^^ and 
described by Holland and Marder^^ who observed it in 
molecular dynamics simulations of silica samples. Exper- 
iments also show that it is possible to obtain atomically 
ffat fracture surfaces in real crystalline materials by using 
sufficiently small and homogeneous loadings. 

The majority of the experiments on cracks produce 
non-planar, rough and branched cracks. The departure 
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from the planar geometry forced by the lattice trapping 
is due to two different contributions. First, the magni- 
tude of the applied loading which influences the energy 
delivered to the crack tip: the higher the value, the larger 
the possibility for the crack to open out of plane branches 
due to the increasing transversal stresses. Second, disor- 
der in the material (equivalently disorder in the breakage 
rule) can drive the crack on non-planar paths. 

Disorder is naturally present in all materials and comes 
from a variety of different sources. Atomic vacancies, in- 
clusions, dislocations and grain boundaries are all sources 
of disorder able to influence the macroscopic response. 
Disorder strongly reduces the effects of anisotropy by 
increasing the probability of deviations from planarity. 
This is the main reason why a high level of disorder was 
included in the simulations presented in this article. The 
easiest way to introduce it (and the way followed by us) 
is to introduce a locally variable fracture energy 7(x) 
which varies according to some well defined distribution: 
although very simple, the uniform distribution accom- 
plishes this task extremely well. If 70 is the fracture 
energy in the absence of disorder, local fracture energies 
can be extracted by a uniform distribution centered on 
7o and ranging between and 270; this ensures that the 
mean fracture energy corresponds to the fracture energy 
in the absence of disorder. 

In the presence of disorder even a static planar test 
crack has different tetrahedra along the crack tip be- 
coming breakable at different loadings, whereas without 
disorder they all became breakable at the same loading 
27o/£^(0). With disorder we defined the reference load- 
ing (corresponding to e = 1) to be the value of at 
which 50% of crack tip tetrahedra in a planar static test 
crack are not breakable. 

In practice modestly lower loadings (e.g. e = 0.7) can 
still lead to crack propagation, as breakage of vulnerable 
tetrahedra along the crack edge leads to stress concen- 
tration at and around more resistant tetrahedra. Higher 
values of e lead to more heavily damaged samples. 



IV. FREE-RUNNING FRACTURE 
SIMULATIONS IN THREE DIMENSIONS 

We now focus on cracks when the constraints of fixed 
crack speed and shape are released. A fixed displacement 
is applied to both top and bottom faces, and a starting 
notch is prepared on the front face. The starting notch 
is long enough to start the simulations in the long crack 
limit. Periodic boundary conditions are applied to the 
side faces. The front and back faces are left stress free. 

Each tetrahedron in the sample was given a pre- 
assigned fracture energy drawn from a broad uniform 
distribution centered on 70 as described in the previous 
section. At each timestep the energy of each tetrahe- 
dron is evaluated, and those in which this exceeds their 
fracture energy are broken by setting the two elastic con- 
stants A = /i = 0. The presence of disorder leads to 



(a) Planar cracks with disorder 
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(b) Non-planar cracks: 300 X 60 X 60 
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(c) Non-planar cracks: 500 X 120 X 
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TABLE I: Measurements of the crack speed v/vr for dif- 
ferent driving energies, in different types of simulations: (a) 
cracks forced to be planar, with the addition of disorder in 
the threshold energy; (b) free cracks with disorder for sam- 
ples 300 X 60 X 60 tetrahedra wide; (c) free cracks with disorder 
for samples 500 x 120 x 120 tetrahedra wide. In all tables, 
the first column is the driving factor e, the second and third 
columns are the crack speed in units of the transverse wave 
speed as measured by the histogram method (second column) 
and the average tip position method (third column). For all 
simulations, vr ~ 0.933 vt. 



breakage of isolated tetrahedra as soon as the simulation 
starts. A condition of connected fracture is imposed by 
allowing only neighbours of already broken tetrahedra to 
break. This simplifies the track of the crack tip in order 
to measure its speed. The condition of connected frac- 
ture prevents the formation of precracks in front of the 
crack tip. In practice, if this condition is released there 
is essentially no difference in the results for moderated 
loadings. We will see in section IVTI that for high loadings 
this can lead to a different morphology for the averall 
fracture process. 

Each simulation starts from the sample relaxed to its 
configuration of minimum energy, and continues until the 
sample is broken into two halves. We wait until the num- 
ber of broken tetrahedra per timestep is reduced to a neg- 
ligible fraction of its peak and then the sample is divided 
into two halves according to the sign of the vertical dis- 
placement of the sites. In some cases the sample remains 
connected by a few isolated tetrahedra, but as these will 
be in a state of anomalously high strain, the displacement 
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FIG. 2: The histogram method for measuring the crack speed. 
The histogram represents the number of broken tetrahedra 
N(l) at distance / from the front side of the sample. As time 
advances, the histogram grows: we can build a set of his- 
tograms corresponding to a set of intermediate positions of 
the crack. The approximate position of the crack front is 
then obtained as the intercept with the dashed threshold. 



of their sites will reflect the displacement of the portion 
of sample to which they are attached. 

That the surface found corresponds to the fracture sur- 
face has been checked using a different method based on 
percolation. The sample can be described as a collec- 
tion of boxes (the tetrahedra) connected on a cubic lat- 
tice. If we choose one tetrahedron on the bottom face we 
can imagine injecting it with some coloured liquid: the 
liquid will spread within the sample and reach the top 




v/v^ 



FIG. 3: Comparison between the results for free-running pla- 
nar cracks (filled circles) and the efficiency description (con- 
tinuous line) . Higher values of the driving energy e correspond 
to lower values of E{v) as described by eq. (|2.3|) . The figure 
shows how this also corresponds to higher crack speeds. The 
dashed line refers to the special case of e = 1. 



face. Broken tetrahedra can act as blocking boxes for 
the liquid, so that the liquid will reach the top face un- 
less a complete fracture surface dividing the sample into 
two halves is retrieved. Complete failure was assured by 
slowly increasing the applied displacement until percola- 
tion between top and bottom was lost, requiring signifl- 
cant extra simulation time. Comparing the two methods 
for a set of 10 simulations a difference of less then 0.05% 
for the set of broken bonds with respect to the whole 
sample was found, and similarly a difference in the flnal 
detected surface of less than 0.9%. 



A. Measurement and selection of Crack Speed 

A flrst series of simulations was built using samples 
with sides 300 x 60 x 60 tetrahedra, with a starting notch 
100 tetrahedra long. In these simulations, the condition 
of flxed crack speed was released and the Griffith criterion 
was used, but the condition of planar ity was maintained 
so that comparison could be made with the efficiency de- 
scription. A sample of three simulations were performed 
for each value of e: above a threshold value of 0.7, all 
simulations led to a completely broken sample. Below 
this value, we found that cracks did not reach the end of 
the sample. 

The constraint of planarity was lifted in the second and 
third series of simulations. The second had parameters 
matched to those of the first, except that larger sam- 
ples of 10 simulations were used to counter the greater 
variability in results. The third series had larger size 
simulations of 500 x 120 x 120 tetrahedra and due to 
computational cost was limited to samples of three sim- 
ulations. 

Crack speeds were measured with two different meth- 
ods. In the histograms method a set of histograms was 
build, reporting the number of broken tetrahedra at each 
distance from the crack notch (see figure (21). Then, the 
intercepts with about 1 /3 of the average value of broken 
tetrahedra per unit of crack advance (dashed line in fig- 
ure) were taken as a measure of the crack front position. 
The value of 1/3 was chosen because that is where the 
histograms are typically steepest. Plotting the position 
against time, the crack speed was retrieved. 

In the average crack tip position method^ the average 
position of the newly-broken tetrahedra in each time in- 
terval considered was retrieved as a function of time, and 
the crack speed was measured from its slope. 

In the first set of simulations, corresponding to planar 
cracks, the values of the crack speed found from both 
methods are fully compatible with the results found in 
our previous work^^ as shown from table ^a) and figure 
IHl The results found are interesting not only for their 
agreement with the results from the efficiency descrip- 
tion, but also because the measured values of v/vr are 
compatible with the crack speeds measured in anisotropic 
materials f^&i in particular with the value of 0.9 vr mea- 
sured for the sample of polymethilmetacrylate prepared 
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FIG. 4: (Color online) (Top) Crack in a sample 300 x 60 x 60 tetrahedra wide (projection), for e = 1. Shading (and color online) 
refers to the depth in the third dimension. The crack appears "fat", thicker than one layer of tetrahedra. A single-layer wide 
section of the sample (bottom) shows that the crack itself is driven slightly off planar. Microbranching is visible, suggesting 
attempted branching as the mechanism which leads to thicker cracks and slows the crack speed. 



with a weak interface/ That value was supposed to be 
the best proof that the limiting crack speed for cracks is 
the Rayleigh speed. The microscopic structure of poly- 
methilmetacrylate is completely different from that of our 
model, however due to the fact that the drop of the ef- 
ficiency at high speeds is dominated by the resonance at 
the Rayleigh speed (which we expect it to be a common 
feature to all materials), we suggest that the limiting 
crack speed for a non-branching crack is the Rayleigh 
speed only in the limit of infinite loading, or e ^ oo. For 
finite loadings, the value of the crack speed measured 
corresponds to that given from the efficiency description. 

The picture is notably different when the condition 
of planarity is released in the second and third series 
of simulations. The disorder dominates over the lattice 
anisotropy, with particularly resistant (or particularly 
weak) tetrahedra prompting the crack to deflect out of 
plane. For low driving forces, the resulting cracks are 
"fat" (figure 21) and the crack speed is considerably lower 
than expected from the efficiency description. Crack 
speeds measured by both methods are reported in table 




FIG. 5: (Color onhne) (Top) Crack in a sample 300 x 60 x 60 
tetrahedra wide (projection), e = 5: branches reach the 
sample's boundaries. (Bottom) Crack in a larger sample, 
500 X 120 X 120 tetrahedra wide (projection), for e = 5. In this 
case branches do not reach the sample boundaries. Shading 
(and color online) refers to the depth in the third dimension. 



IJb). Low values of e give a crack speed which is lower 
than expected from the efficiency description, but com- 
patible with some of the measurements obtained both in 
simulations and in experiments as explained in section 
n For higher values of e, the results of the two types of 
measurements for the crack speed diverge, which we at- 
tribute to the influence of crack branching. The average 
of the crack tip position is an average of the new bro- 
ken tetrahedra, which includes tetrahedra which break 
along branches behind the crack front. As a result, the 
average is affected by a systematic error which lowers 
the value from that of the crack front. Where the speed 
measurements differ significantly, we take the histogram 
method to give the true speed of the crack front and the 
divergence of the methods to be an indicator of crack 
branching. 

For the highest values of e in the set of smaller simu- 
lations, the histogram method gives a value of the crack 
speed above that predicted from the efficiency function 
and, for e = 5, above the Rayleigh speed. This unex- 
pected result is a consequence of the high level of dam- 
age of the sample due to the high level of the driving 
force. Direct inspection (see fig. 0) shows that the sam- 
ple is broken, with branches reaching the sample's top 
and bottom boundaries. The anomaly disappears upon 
increasing the size of the sample. Simulations for sam- 
ples 500 X 120 X 120 tetrahedra wide show, for low values 
of e, results compatible with those of the smaller samples 
and, for all values of e, crack speeds lower than the values 
given by the efficiency function (see table ^c)). Direct 
inspection (see fig. 0) shows that although the damage 
is still heavy, in this case branches no longer reach the 
sample boundaries. 

How do these results compare with those of other sim- 
ulations and experiments? Experiments are usually per- 
formed with carefully controlled loading. The load is 
slowly increased and put just above the level beyond 
which the macroscopic fracture develops. This cor re- 
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FIG. 6: Sequence of frames for a typical crack growing in a sample 500 x 120 x 120 tetrahedra wide. The gray shading refers 
to the height with respect to the crack notch plane. In this case e = 3. 



spends to a low level of e. This can help explain the 
compatibility of these values of the crack speed with 
those measured in experiments. With respect to the re- 
sults of the planar case however, crack speeds are lower. 
Although no major branches develop as testified by the 
agreement in the values of the crack speed for the two 
type of measurement, such drop in crack speed with re- 
spect to the efficiency prediction must be connected with 
the mechanism of attempted branching which is respon- 
sible for the "fat" appearence of these cracks. For e = 1, 
the average width of a crack is 2.35 tetrahedra, indicat- 
ing that more energy is needed for the crack to advance 
than expected from the case of a the fully planar crack. 

Higher values of e do not correspond to the experimen- 
tal set up for the measurement of crack speeds. For such 
values, the crack speed approaches the values given by 
the efficiency description. The growth of the crack speed 
for increasing loading is similar to that measured in other 
simulations^^'^^. The high driving force regimes are those 
which give rise to a non-zero roughness exponent as we 
will see below. 



B. The shape of advancing cracks 

Figure shows a sequence of frames of an advancing 
crack for e = 3, up to its final state, when the sample is 
fully broken. Due to the magnitude of e, a considerable 
level of branching is visible. Straight cracks are obtained 
for values of e < 1 (see later fig. [T^ . 

At first sight, a reasonable description for the phe- 
nomenon is that of a main crack from which minor 
branches develop during the dynamics. This is not a com- 
plete picture. Side views of the same final fracture (see 
fig.Q) reveal that the crack is made by a set of connected 



branches. The crack tip splits into two or more branches 
which try to avoid each other, and force the crack to 
advance in a non linear fashion. The particular fracture 
shown is also characterized by two main branches running 
along much of the sample, and comparing the different 
panels of figure [7| it can be seen that both of these are 
part of the final fracture surface. Most branches die out 
as they head towards the sample's boundaries. Those 
branches that mantain their distance from the sample's 
boundaries manage to travel through the sample build- 
ing up its backbone. The development of branches and 
then of the backbone, is clearly controlled by the sam- 
ple's boundaries which act as a guide to the branching 
process and lead the whole crack in the forward direction. 



V. ROUGHNESS OF FRACTURE SURFACES 

An example of a final fracture surface is shown in figure 
IHl The surface does not appear flat, and its roughness 
can be quantified through the roughness exponent. For 
two points separated along the direction of (global) crack 
proagation we have 



{h{x) - h{x')Y 



(X \x ■ 



and similarly we define ("^ for the scaling of height fluc- 
tuations along the direction of the (global) crack edge. 
These scaling laws can equivalently be probed by spatial 
power spectra, so that for a cut along the x-direction we 
expect 



h{kxf) (X \kx 



-i-2Cx 



The scaling is expected to apply from local lengthscales 
up to of order the (smallest) dimension of the sample. 
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We find the measured roughness exponent varies sys- 
tematically with the strength of crack driving e, whilst 
there is relatively little difference between measurements 
of the exponent in different directions or by different 
methods. To limit the influence of the boundaries and 
of the starting notch, we analysed a region limited to the 
central 80 x 60 tetrahedra for the 300 x 60 x 60 samples, 
and 120 x 120 tetrahedra for the 500 x 120 x 120 samples, 
equidistant between the final boundary and the end of 
the starting notch. 

Correlation functions for the 300 x 60 x 60 samples, 
shown in figure [Ul appear ordered in e for their slope. 
For increasing values of the driving energy, a region of 
constant slope close to the origin develops in the cor- 
responding correlation function. The extension of this 
region is larger, the higher the value of e. From its slope 
we have measured the roughness exponent corresponding 
to each value of the driving force. The roughness of the 
surface grows up to a limiting value, as shown by the 
common slope of the highest curves in both figures. The 




FIG. 7: (Color online) Slices of the fracture of figure El viewed 
from the side. Each slide corresponds to 1/8 of the sample. 
Shadings (and colors online) refer to the depth in the third 
dimension. 








FIG. 8: Final fracture surface corresponding to a driving en- 
ergy e = 3. 




1 5 20 




1 5 



FIG. 9: Log- log plots of the height-height correlation func- 
tions along the x direction (on the top), and the z direction 
(on the bottom) for the 300 x 60 x 60 samples. For comparison, 
also one of the correlation functions for the 500 x 120 x 120 
samples is shown. The above slope corresponds to roughness 
exponent C = 0.45. 

loss of slope of the highest curve corresponding to e = 5 
can be explained due to the approach of branches to the 
sample's boundaries as described in the previous section. 

Figure irni shows the power spectra of the same cuts on 
the surface. Although the common roughness exponent 
is retrieved, the characteristics of how this common slope 
builds up are less clear, and data are more scattered. 

Results for the roughness exponent for both sets of 
samples are reported in figures ITTl and IT^ Results from 
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FIG. 10: Log-log plots of the height power spectra for cuts 
along the x direction (on the top), and the z direction (on 
the bottom) for the 300 x 60 x 60 samples. The above slope 
corresponds to ^ = 0.45. 



the power spectra are compatible with the real space 
ones, although again more scattered. These results show 
that the roughness exponent increases quite sharply from 
a value close to zero, which we will see below can be con- 
nected with a logarithmic scaling, up to a value between 
0.4 and 0.5 for both the x and z directions. Compari- 
son between the values of the roughness exponent along 
the two space directions shows a slight difference which 
might be attributable to boundary effects, although data 
are insufficient to draw a definite conclusion. 

For the lowest values of e at which we could propagate 
cracks, our roughness data are better described by a log- 
arithmic roughness law as shown by figure ^1 This is 
interesting because it matches calculations by Ball and 
Larralde^^ and subsequent calculations by Ramanathan, 
Ertas and Fisher^^ based on continuum elastic fracture 
mechanics for cracks at the threshold of propagation, 
as well as supporting experiments reported in Ref. Q. 
In continuum elastic fracture mechanics these cracks are 
quasi-static, which we know from the efficiency descrip- 
tion is not the case for a structured material, whilst in 
the supporting experimentsii the overall propagation was 
kept slow but locally crack acceleration could (and did) 
occur. Thus it appears that the logarithmic law may ap- 
ply rather generally to cracks propagating marginally in 
three dimensions, without restriction to zero speed. 

The increase in the fracture roughness can be visually 
connected to the branching process by comparing these 
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FIG. 11: Change of the roughness exponent from the corre- 
lation functions, with the driving factor e for cuts along the 
X direction (on the top) and the z direction (on the bottom). 
Results for the 300 x 60 x 60 samples are in black, those for 
the 500 X 120 x 120 samples are in light gray. 



results to the sequence of figure The sequence shows 
that cracks are flat for low values of e within the limits of 
what disorder allows. Macroscopic branches appear at a 
value of e around 1.4 — 1.6 which corresponds to the in- 
crease of the roughness exponent towards its limit value. 
The result should be compared with what observed by 
Sharon, Gross and Fineberg^^ in their experiment on the 
branching instability. The shape of branches appeared to 
be compatible with a power law shape with an exponent 
of about 0.7. Our results suggest that the appearence 
of branching increases the roughness to a value which in 
our case is between 0.4 and 0.5, which corresponds to the 
one measured at short length scales. 



VI. SIMULATION OF DISCONNECTED 
FRACTURES 

The limit value for the roughness exponent at about 
0.45 corresponds to the one measured in molecular dy- 
namics simulations by Nakano, Kalia and Vashishta^^i 
in microcrack advance, where they suggest microcrack 
coalescence as the mechanism leading to higher values 
~ 0.75.^^ As all our simulations above have been per- 
formed imposing the connection of the advancing crack, 
which prevents the formation of multiple cracks within 
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Height-height correlation function 


Power spectrum 


e 


c. 




c. 




1.0 


0.153 ±0.008 


0.265 ±0.006 


0.058 ± 0.033 


0.030 ± 0.060 


2.0 


0.330 =b 0.002 


0.419 ±0.003 


0.280 ± 0.050 


0.402 ± 0.087 


3.0 


0.339 ± 0.003 


0.445 ± 0.005 


0.341 ± 0.042 


0.661 ± 0.050 


4.0 


0.473 zb 0.005 


0.387 ±0.001 


0.525 ±0.011 


0.398 ± 0.054 


5.0 


0.418 ± 0.002 


0.452 ± 0.004 


0.509 ± 0.043 


0.526 ± 0.094 



TABLE II: Roughness exponent for different values of e in the case of non- connected fractures, as discussed in Section IVll All 
these results have been obtained from simulations of samples 500 x 120 x 120 tetrahedra wide. The second and third column 
correspond to the roughness exponent measured from the height-height correlation functions along cuts in the x and the z 
directions. The fourth and fifth columns correspond to the roughness exponent measured from the power spectra. 
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FIG. 12: Change of the roughness exponent from the power 
spectra, with the driving factor e for cuts along the x direction 
(on the top) and the z direction (on the bottom). Results for 
the 300 X 60 X 60 samples are in black, those for the 500 x 
120 X 120 samples are in light gray. 



the sample, we could expect to have the same coexis- 
tence of the two values for the roughness exponent if 
we release this additional condition. Its removal allows 
the creation of diffuse damage (commonly referred to as 
"dust"), composed of isolate broken tetrahedra scattered 
throughout the sample. The level of such dust increases 
with e, as the number of breakable tetrahedra increases 
with e. In fig. the case for e = 5 is shown. The 
outcome has been cleaned of the dust which, due to the 
extreme high value of e, sums up to a number of broken 
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FIG. 13: Linear-log plots of the height-height correlation 
functions along the x direction (on the top), and the z di- 
rection (on the bottom) for the 300 x 60 x 60 samples. In 
both plots, the lowest curves corresponding to the lowest e 
are straight, corresponding to a logarithmic scaling. 



bonds equivalent to those shown, but scattered through 
the whole sample. The final fracture is characterized by 
large jumps due to microcracks which have connected 
during the dynamics, suggesting a similar interpretation 
to that of Nakano et al. Measurements of the correspond- 
ing roughness exponent for the 500 x 120 x 120 samples 
however, reported in tabled do not show any relevant 
difference from those of connected fractures. This may 
just reflect the length scale at which microcrack coales- 
cence develops in our simulations being larger than the 
length scale over which our roughness measurements are 
performed. To check any increase in the roughness expo- 
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FIG. 14: Shape of cracks in a sample 300 x 60 x 60 tetrahedra 
wide (projections) for increasing values of e. Values shown are 
from the top: e = 0.8, 1.0, 1.4, 1.6, 2.0, 3.0, 4.0, 5.0. Shading 
refers to the depth in the third dimension. 



nent associated with the lengthscale of crack coalescence 
evident in fig. ^| an increase in the system size would 
certainly be necessary. 



VII. CONCLUSIONS 

The release of the condition of fixed crack speed has 
been the first test for the efficiency description of the 
selection of the crack speed. The efficiency description 




FIG. 15: Sideview of the final fracture in the case of non- 
connected fracture as discussed in Section IVII for a 500 x 
120 X 120 tetrahedra wide sample and e = 5. Shading refers 
to the depth in the third dimension. 



shows, in contrast to the theoretical description of the 
continuum theory, that the Rayleigh speed is the ultimate 
crack speed only in presence of an infinite loading. For fi- 
nite loading applied, the crack speed selected follows the 
efficiency description in presence of consistent anisotropy 
or lattice trapping. That the terminal crack speed would 
be the Rayleigh speed only for an infinite loading had al- 
ready been suggested by Xu and Needleman.^^ However, 
here this result comes out naturally from the shape of 
the efficiency, due to the Griffith criterion (|2.2|) . 

When the fracture geometry is released from planar, 
the terminal crack speed in our simulatiosn is much lower 
than expected from the efficiency description, and corre- 
sponds to the range of crack speed measured in experi- 
ments. The analysis of the crack shape has shown that 
even at the lowest loadings, cracks attempt to branch 
and this tends to thicken the crack. Hence, the amount 
of energy which is delivered into fracture work per unit of 
crack advance has to increase, and the energy radiated 
has to decrease. This inevitably leads to a change in 
the maximum speed achieved. Although the basic idea is 
clear, more work is required to relate the effective amount 
of energy radiated with the branching mechanism and 
hence with a criterion of speed selection for non-planar 
cracks. 

When the energy is sufficient, branching becomes a 
macroscopic phenomenon and appears to be the basic 
mechanism through which the roughness exponent of a 
crack surface builds up. The dynamics of the macroscopic 
branches determine the final fracture shape: branches 
build the backbone of the whole crack. During the 
crack advance branches try to avoid each other whilst 
the boundaries of the sample lead them into the forward 
direction, creating a preferential direction for their ad- 
vance. All the branches that greatly deviate from this 
direction die out and do not belong to the final fracture 
surface. 

When the driving energy is low, the self-affine prop- 
erties of the crack surface are also compatible with the 
scenario of logarithmic roughness suggested in Ref. 
and El The condition of quasi-static cracking is explicit 
in the latter in terms of the fracture energy supplied to 
the crack being barely in excess of the fracture toughness. 
This however does not correspond simply to having a 
crack with little kinetic energy. On the contrary, because 
of the efficiency description, low values of e correspond to 
cracks which barely have energy to advance, but which 
travel at high speed. This could clarify how a logarithmic 
roughness predicted for a quasi-static crack can be found 
in some of our simulations of dynamic cracks. 

The roughness exponent grows rapidly with the driving 
energy e, towards a maximum of ~ 0.45 corresponding to 
that measured in experiments on short length scales. Re- 
cent experiments^^'^^ have shown that fracture surfaces 
could be affected by anomalous scaling. This in turn 
has been interpreted^^ as the mark of the anisotropy of 
the roughness exponent measured at large lenght scales. 
Our simulations show a slight anisotropy of the final frac- 
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ture surface, which might be an indication of similar be- 
haviour on short length scales, and calls for further at- 
tention in both simulations and experiments. 

We have not been able to corroborate the suggestion 
of Ref. 42 that the higher roughness exponent for large 
scales results from microcrack coalescence, but this may 
be due to the limited range of lengthscales we have been 
able to access in these simulations. We estimate we would 
require to gain one to two orders of magnitude in distance 



range to resolve this issue with the present model. 
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